Hubbard models for materials with nonlocal Coulomb interactions: graphene, silicene 

and benzene 
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We report on a method to approximate a generahzed Hubbard model with nonlocal Coulomb 
interactions by an effective Hubbard model with on-site interactions U* only. The effective model 
is defined by the Peierls-Feynman-Bogoliubov variational principle. We find that the local part of 
the interaction U is reduced according to = — V, where V is an average of nonlocal Coulomb 
interactions weighted according to derivatives of nonlocal charge correlation functions from the 
effective model. For the examples of graphene, silicene and benzene we show that the nonlocal 
Coulomb interaction can decrease the effective local interaction by more than a factor of 2. 

PACS numbers: 72.80.Rj; 73.20.Hb; 73.61.Wp 



Low dimensional sp-electron systems like graphene 
3], systems of adatoms on semiconductor surfaces, such 
as Si(lll):X with X=C, Si, Sn, Pb Q, Bechgaard salts 
or aromatic molecules [H, [6| feature simultaneously strong 
local and nonlocal Coulomb interactions. In graphene for 
instance, the on-site interactions U/t ^ 3.3 , the near- 
est neighbor Coulomb repulsion F/t ~ 2 as well as fur- 
ther sizable nonlocal Coulomb terms exceed the nearest 
neighbor hopping t = 2.8 eV Considering on-site 
interactions U/t ~ 3.3 alone would put graphene close 
to the boundary of a gapped spin-liquid [7[ , which could 
be even crossed by applying strain on the order of a few 
ercent It is currently unclear, whether ^ or not 
nonlocal Coulomb interaction stabilize the semimetal- 
lic Dirac phase in graphene. To rephrase the problem: 
It is unclear which Hubbard model with strictly local 
interactions would yield the best approximation to the 
ground state of graphene. To judge the stability of the 
Dirac electron phase in graphene but also to understand 
Mott transitions on surfaces like Si:X (111), a quantita- 
tive well defined link from models with local and nonlocal 
Coulomb interactions to those with purely local interac- 
tions is desirable. 

In this letter, we present a method to map a gener- 
alized Hubbard model with nonlocal Coulomb interac- 
tions onto an effective Hubbard model with on-site in- 
teractions U* only. Using the examples of graphene, 
silicene and benzene we show that nonlocal terms can 
significantly reduce the effective on-site interaction. For 
graphene, this reduction is U* /U ^ 0.45. Thus, nonlo- 
cal Coulomb interactions are found to stabilize the Dirac 
electron phase against spin-liquid and antiferromagnetic 
instabilities. We begin with a description of the theoret- 
ical framework and afterwards present its application to 
graphene, silicene and benzene. 



The starting point is the extended Hubbard model 
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where tij are the hopping matrix elements and U and 
Vij are the local and nonlocal Coulomb matrix elements, 
respectively. The goal is to map the Hamiltonian ([1]) onto 
the effective model 
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The effective on-site interaction /7* shall be chosen such 
that the canonical density operator p* = 1/Z*e~^^ of 
the auxiliary system, where Z* = Tr{e~^^ } is the 
partition function, approximates the exact density op- 
erator p derived from H as close as possible. This re- 
quirement leads to the Peierls-Feynman-Bogoliubov vari- 
ational principle [lol-llij for the functional 
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where ^* = — InZ* is the free energy of the auxiliary 
system. (...)* denotes thermodynamic expectation val- 
ues with respect to the auxiliary system: {H — i^*)* = 
Trp*(i^-i^*). In the case of p* = p the functional <l[p*] 
becomes minimal and coincides with the free energy. The 
optimal [/* is thus obtained for minimal = 4> [[/*]: 



S^*<|[/7*] =0. 
By evaluating Eq. (|4]) one finds 

u* = u + lJ2Vi 
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Figure 1. (Color online) Illustration of physical process un- 
derlying Eq. (|5]). Wavy lines illustrate Coulomb interactions, 
(a) An electron in the extended Hubbard model hopping from 
a doubly occupied place to an empty place, gaining an energy 
U — V. (b) The same situation in the effective model leads to 
an energy gain of [/*. 



This rule presents a central result of this letter and has an 
intuitive physical interpretation (Fig. [T]): Increasing the 
on-site term /7* reduces the double occupancy (n^^n^i)* 
and pushes away electrons approaching an already occu- 
pied site i = to neighboring sites. In case of purely local 
Coulomb interactions there is a Coulomb energy gain of 
[/* upon suppressing the double occupancy (Fig. (Tb). 
When there are, however, nonlocal Coulomb interactions 
with surrounding lattice sites j, the displaced electrons 
raise the energy of the system by terms proportional to 
Voj. For a simple case of two electrons on one site this 
process is depicted in Fig. [TJi. In this case, it is obvi- 
ous that the Coulomb energy gain due to the electron 
displacement in the full and the auxiliary model become 
energetically equivalent for /7* = /7 — V. In general, this 
energy gain depends both on to which sites the charge 
density is displaced due to the local Coulomb interaction 
and on how strong the nonlocal Coulomb terms are. 

For a translationally invariant system, the local part 
of the interaction U is reduced according to /7* = /7 — 
where 



Oj 
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9f/*(notno;)* = -X;^yo,a^t/*(^ot^ja)*. Thus, V is a 
weighted average of the nonlocal Coulomb interactions. 
It is reasonable that V is positive (repulsive) in most sit- 
uations that correspond to real materials. Then the non- 
local Coulomb interaction reduces the effective on-site 
interaction and therefore stabilizes the Fermi sea against 
transitions e.g. to a Mott insulator. Nevertheless, models 
with negative V can easily be constructed and it remains 
to be seen which real materials lead to negative (attrac- 
tive) V. 

Under the assumption that an increasing /7* displaces 
electrons only to next neighbors, we find du* (notno;)* = 
—Nndu* X]cr(^ot^icr')*5 where Nn is the coordination 
number. Equation. (|5]) then yields 



/7* = /7-K)i. 



(7) 



The conservation of the total electron number N 
leads to the sum rules Xl7cr(^ot^jcr)* = N/2 and 



This result is appealingly simple and gives an estimate 
for the effective Coulomb interaction, without the need 
of numerical calculations. It follows, however, from a 
severe approximation. The following numerical calcula- 
tions show that Eq. ([7j) nevertheless leads to values close 
to the exact ones (shown in table H]). 

When the approximation that electrons are only dis- 
placed to next neighbors is dropped, the derivatives of 
the correlation functions have to be calculated explicitly. 
This can be done approximately within the dynamical 
mean field theory [13| and diagrammatic extensions like 
the Dual-Fermion approach [14]. In certain cases also nu- 
merically exact calculations of the nonlocal charge cor- 
relation functions for instance by means of determinant 
quantum Monte Carlo (DQMC (HI) or density-matrix 
renormalization group methods (e.g. (HI) are possible. 

In the following, we consider graphene and benzene, 
where the computational effort of the DQMC method is 
manageable for sufficiently small lattice sizes. Therefore, 
we used the DQMC implementation "QUantum Electron 
Simulation Toolbox" (quest 1.3.0 [13) on a super cell 
to obtain the charge correlation functions that enter Eq. 
(|5j). For small systems a different implementation [18| 
was used to verify the results of the QUEST package. The 
Hubbard model with less than 8-9 sites can also be solved 
by exact diagonalization (ED). In this case a comparison 
with data obtained with DQMC shows excellent agree- 
ment. To overcome the statistical noise from the DQMC 
calculations, we fit the correlation functions with poly- 
nomials of rank 4 in [/* and evaluate the derivative an- 
alytically. To evaluate Eq. (j5j) in a numerically stable 
way for extended systems like graphene, a cutoff Tc is in- 
troduced. All derivatives of the correlation functions for 
which the distance between site i and j is bigger than Vc 
are neglected. A reasonable convergence is found for a 
cutoff at the fourth next neighbor. 

To calculate [/* for realistic systems, we introduce val- 
ues for the Coulomb interactions in the original model 
defined by Eq. ([1]). For graphene and silicene these val- 
ues are calculated with the constrained random phase 
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Figure 2. (Color online) Derivatives of the correlation func- 
tions (no^rija) with respect to U* for graphene (16x16 unit 
cells) calculated with DQMC for inverse temperature 13 = 9t 
and U — 2t. Each circle corresponds to one carbon atom. The 
shaded area depicts the region of nearly vanishing derivatives. 
The cutoff radius is Vc The thick drawn circles indicate the 
lattice site with index i = 0. 

approximation (cRPA) [IqI Hke in [l|. Eor benzene, we 
use values from [2o| , which are obtained by fitting U and 
t to experimental spectra and calculating Vij by Ohno in- 
terpolation [21]. For all systems the values of the initial 
Coulomb interactions are given in table HI 

For a honeycomb lattice at half filling the 
derivatives of the correlation functions du* {riQ^njcr'Y for 
U* = 2t are shown in Fig. [2l In this particular case, 
du* (^ot^jcr')* changes sign with both sublattice and spin 
indices. Generally, (no^nj^)* | with opposite spins 
(right panel) exceeds the equal spin case \du* {no^nj^y\ 
(left panel). The derivatives decrease with the distance 
between the sites i = (thick drawn circles in the mid- 
dle) and j. According to Eq. (|6]), the derivatives of the 
correlation functions are the weight of the average of the 
Voj to calculate /7*. Thus, the oscillating behavior of 
the correlation function with respect to the sublattice as 
well as the spin direction speeds up the convergence of 
Eq. with the cutoff Tc. 

An exemplary result of the evaluation of Eq. (|5j) for 
different rc is shown in Fig. [3l (See [22^ for compu- 
tational details.) The error bars are determined by the 
statistical errors of the DQMC calculations. As could be 
already expected from the du* {no^n jcr')* -terms depicted 
in Fig. O the nearest neighbor Coulomb interaction has 
the strongest impact on the renormalization of the on- 
site term /7*. Thus, the approximation of neglecting all 
non nearest neighbors terms is quite reasonable, here. 
Contributions from the same sublattice as the site i = 
tend to increase [/* and those from the other sublattice 
decrease t/*. It is obviously sufficient to introduce the 
cutoff after the fourth nearest neighbors. 

The resulting values of the effective local Coulomb in- 
teraction /7* for graphene, silicene and benzene are sum- 



Figure 3. (Color online) Effective Coulomb interaction U* for 
graphene against the position of the cutoff used to evaluate 
Eq. {O. 

marized in table HI The local Coulomb interaction is de- 
creased by a factor of larger than two in all cases. For 
both, graphene and silicene the renormalized on-site in- 
teractions /77t = 1.6 ± 0.2 and /7*/t = 2.0 ± 0.3 are 
far away from the transition to a gapped spin liquid at 
/t = 3.5 [7]. The Dirac semimetal phase is thus stabi- 
lized by the nonlocal Coulomb interactions. 

We obtain the strongest renormalization of the on- 
site interaction for benzene. This is mostly due to the 
different ratio between local and nonlocal Coulomb in- 
teractions in benzene, Vqi/U = 0.72, as compared to 
Voi/U = 0.56 for graphene Q or Voi/U = 0.55 for sil- 
icene. 

Table I. First three rows: Coulomb matrix elements obtained 

with cRPA (tgraphene = 2.80 cV, tgilicene = 1-14 cV, tbenzene = 

2.54 eV). Last three rows: Effective local Coulomb matrix 
elements with and without the approximation that electrons 
are only displaced to next neighbors and the factor by which 
the local Coulomb interaction is decreased. 





Graphene 


Silicene 


Benzene 


u/t 


3.63 


4.19 


3.96 


{Voiyo2)/t 


2.03,1.45 


2.31,1.72 


2.83,2.01 


(H3,K)4)/t 


1.32,1.14 


1.55,1.42 


1.80 


uyt 


1.6 ±0.2 


2.0 ±0.3 


1.2 


{U-Voi)/t 


1.6 


1.9 


1.1 


uyu 


0.45 ±0.05 


0.46 ±0.05 


0.3 



Naturally, the question arises how accurate the effec- 
tive model reflects the physical properties of the orig- 
inal model. The phase diagram of the extended Hub- 
bard model on the honeycomb lattice includes an anti- 
ferromagnetic (AF), a semimetal (SM) and a charge den- 
sity wave (CDW) phase [2^, while the Hubbard model 
with strictly local interactions only features the first two 
phases. If the parameters of the extended model are 
clearly inside the AF or the SM phase, the effective model 
will likely approximate the physical properties of the orig- 
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Figure 4. (Color online) Correlation as functions of the screen- 
ing for the extended Hubbard model (full model, continuous 
lines) and the effective Hubbard model (broken lines) for ben- 
zene. The left panel shows the spin correlation function (Sl^) 
and the right panel shows the density correlation function 
{pl^). {Sz^) is virtually the same for the effective and full 
model. The parameters for the full model are U = 10.06 eV, 
tfuii = 2.539 eV and the non local Coulomb interaction V{e) 
is calculated by Eq. dH}. The local interaction for the effective 
model U*{s) is calculated by (|5)), while teff. — tfuii- 



inal model quite well. If the full system is, however, close 
to the CDW phase, a failure of the auxiliary model to 
provide a physically correct description of the original 
situation can be expected. 

We illustrate this expectation with the example of 
modified benzene. In this modified model, the nonlocal 
Coulomb interaction 



U 



a/1 + {aenjy 



(8) 



with a = U/e'^ is tuned by an additional variable screen- 
ing e ranging from to oo. e = oo corresponds to purely 
local interactions and e = ultimately nonlocal inter- 
actions with matrix elements not decaying with distance 
between the involved electrons. A comparison of the spin 
(Sl^) = {{ui^ — nii){nj^ — Uji)) and the density correla- 
tion functions (p*-^) = {{ui^ -\-nii){nj^ + ^^e ex- 
tended and the auxiliary local Hubbard model are shown 
in Fig. HI The correlation functions shown here have 
been calculated by exact diagonalization for, both, the 
full and the effective model. For e = and e ^ oo 
(non-interacting limit) the correlation functions of the 
effective and full model coincide as they should. CDW 
physics would manifest in (p*-^ ) and here we find indeed 
some differences of (p*-^) for the effective and the auxil- 
iary model for intermediate screening (e ~ 1). However, 
nearly no deviation of (Sl^) between the full and effec- 
tive model is found. We thus expect that transitions into 
phases like an AF insulator (or a Mott insulator) will be 
very well described by the effective model. 

In conclusion, a systematic map from a generalized 
Hubbard model with nonlocal Coulomb interactions to 



an effective Hubbard model with strictly local Coulomb 
interactions U* is derived. The physical properties of the 
effective model reflect the original system nicely, espe- 
cially regarding spin related properties. We find that the 
nonlocal Coulomb interactions can significantly renor- 
malize the effective on-site interaction /7* as compared 
to the original local U. In the cases of graphene and sil- 
icene our calculations yield U* /U < 0.5. Thus, the nonlo- 
cal Coulomb interactions stabilize the Dirac semimetallic 
phases in these materials against transitions to a gapped 
spin liquid or an antiferromagnetic insulator. In defec- 
tive graphene or at edges local Coulomb interactions can 
lead to the formation of magnetic moments [sl, [25|, [26| . 
When describing these situations in terms of the Hubbard 
model, the value of /7* = IM obtained here should be 
used. For all materials studied, here, ^ U — Vbi, Eq. 
([7]), yields a qualitatively correct estimate of the effective 
local interactions without the need for numerical calcu- 
lations. We therefore speculate that nonlocal Coulomb 
interactions will significantly weaken local correlation ef- 
fects in sp-electron materials, in general. 
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